home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Languguage OS 2
/
Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO
/
language
/
parallax
/
more_exa.tar
/
more
/
Graphics
/
raytrace.p
< prev
next >
Wrap
Text File
|
1991-11-18
|
59KB
|
1,732 lines
(***************************************************************************)
(* *)
(* Pixelparallel Raytracing-Algorithm *)
(* *)
(* *)
(* Author: Sabine Liebelt *)
(* File: trace.p *)
(* Language: Parallaxis *)
(* *)
(***************************************************************************)
SYSTEM sys_trace;
CONST
MaxSph = 100; (* maximal number spheres *)
MaxPoly = 100; (* maximal number polygone *)
MaxVertices = 4; (* maximal number polygon-corner *)
MaxLights = 5; (* maximal number Lights *)
MaxDepth = 3; (* maximal depth of raytracing *)
MinWeight = 0.001; (* minimal weighting factor *)
Size = 300; (* image-size *)
RayEps = 0.0001; (* factor of variation *)
TYPE
RGB = ARRAY [1..3] OF REAL; (* values between 0 and 1 *)
Vec = ARRAY [1..3] OF REAL;
Ray = RECORD
p, d: Vec;
END;
Sphere = RECORD
center: Vec;
radius: REAL;
kdr: REAL; (* coefficients for Shading *)
ksr: REAL;
shine: REAL;
kst: REAL;
eta: REAL;
color: RGB;
END;
Polygon = RECORD
normal: Vec;
d: REAL;
xy: BOOLEAN; (* TRUE: projection on xy-plain *)
(* usefull *)
xz: BOOLEAN; (* TRUE: projektion on xz-plain *)
(* usefull *)
side: INTEGER; (* 1:normalvector has pos. z-component*)
(* -1:normalvector has neg. z-component*)
vertices: ARRAY [1..MaxVertices] OF Vec;
vcount: INTEGER;
kdr: REAL; (* coefficients for Shading *)
ksr: REAL;
shine: REAL;
kst: REAL;
eta: REAL;
color: RGB;
END;
SphList = ARRAY [1..MaxSph] OF Sphere;
PolyList = ARRAY [1..MaxPoly] OF Polygon;
ObjTyp = (none, sph, poly);
Object = RECORD
typ: ObjTyp;
poly: Polygon;
sph: Sphere;
END;
Viewer = RECORD
pos: Vec; (* point of view *)
at: Vec; (* direction of view *)
up: Vec; (* direction above *)
angle: REAL; (* view-angle in degree *)
END;
Light = RECORD
pos: Vec;
intensity: RGB;
END;
LightList = ARRAY [1..MaxLights] OF Light;
string = ARRAY [1..50] OF CHAR;
CONFIGURATION Field [1..Size],[1..Size];
CONNECTION;
SCALAR
eye: Viewer; (* viewer *)
depth: INTEGER; (* depth of recursion *)
sphcount: INTEGER; (* spherecounter *)
polycount: INTEGER; (* polygoncounter *)
sl: SphList; (* spherelist *)
pl: PolyList; (* polygonlist *)
lightcount: INTEGER; (* lightcounter *)
ll: LightList; (* lightlist *)
amb: RGB;
background: RGB;
picfile, objfile, inputf: string;
i, j: INTEGER; (* counter *)
pixfield: ARRAY [1..Size],[1..Size] OF RGB;
m: REAL; (* max. component of current *)
(* colorvector *)
aktu_color: RGB; (* color of current pixel *)
VECTOR
iray: Ray;
color: RGB;
Vecpoly: Polygon;
PROCEDURE strcat( SCALAR first , second : string ) : SCALAR string ;
(***************************************************************************)
(* concatenates two strings *)
(***************************************************************************)
SCALAR
i , j : INTEGER;
BEGIN
i := 1 ; j := 1 ;
WHILE first[i] <> CHR(0) DO INC(i) ; END ;
WHILE second[j] <> CHR(0) DO first[i] := second[j] ; INC(j) ; INC(i) ; END ;
first[i] := CHR(0) ;
RETURN first ;
END strcat ;
(********* READ DATA AND CREATE DATASTRUCTURE ******************************)
PROCEDURE createscene ();
(***************************************************************************)
(* *)
(* subprocedures: *)
(* get_viewer, get_light, get_background, get_poly, get_patch, get_sphere,*)
(* get_material, Read_comment *)
(* *)
(* global variables: *)
(* *)
(* function: *)
(* read viewer, object and light-data. *)
(* If no viewer exist, then program will be terminated. If more than one *)
(* viewer exist, then the last one would be used. *)
(* *)
(* Inputfile has to be in NFF-format : *)
(* *)
(* # Viewpoint *)
(* v from pos_x pos_y pos_z *)
(* at at_x at_y at_z *)
(* up x y z *)
(* angle alpha *)
(* hither dist *)
(* resolution x y *)
(* # Lights *)
(* l pos_x pos_y pos_z r g b *)
(* # Backgroundcolor *)
(* b r g b *)
(* # Material *)
(* r g b kdr ksr shine kst eta *)
(* # Spheres *)
(* K cen_x cen_y cen_z rad *)
(* # Polygones *)
(* P num.Corners x y z x y z ... *)
(* *)
(***************************************************************************)
SCALAR
view: BOOLEAN; (* TRUE = viewer known *)
(* FALSE = viewer still unknown *)
back: BOOLEAN; (* TRUE = background known *)
(* FALSE = background still unknown *)
material: BOOLEAN; (* TRUE = surface known *)
(* FALSE = surface still unknown *)
mat: RECORD (* current material *)
color: RGB;
kdr: REAL;
ksr: REAL;
shine: REAL;
kst: REAL;
eta: REAL;
END;
t: string; (* s = sphere *)
(* l = light *)
(* p = polygon *)
(* v = viewer *)
PROCEDURE get_viewer();
(***************************************************************************)
(* *)
(* calling procedure: *)
(* createscene *)
(* *)
(* global variable: *)
(* eye *)
(* *)
(* function: *)
(* to read viewerdata from file *)
(* *)
(***************************************************************************)
SCALAR
i: INTEGER; (* counter *)
next: string;
resx, resy: INTEGER;
hither: REAL;
BEGIN
IF view THEN
WriteString ("WARNING: more than one viewer specified.");
WriteLn;
WriteString (" use last one.");
WriteLn;
ELSE
view := TRUE;
END;
IF Done THEN
ReadString (next);
END;
IF STRCMP (next, "from") = 0 THEN
(* read position *)
i := 1;
WHILE Done & (i <= 3) DO
ReadReal (eye.pos[i]);
INC (i);
END (* while *);
ELSE
WriteString ("ERROR: viewerposition not correct.");
HALT;
END;
IF Done THEN
ReadString (next);
END;
IF STRCMP (next, "at") = 0 THEN
(* read direction of view *)
i := 1;
WHILE Done & (i <= 3) DO
ReadReal (eye.at[i]);
INC (i);
END (* while *);
ELSE
WriteString ("ERROR: direction of view not correct.");
HALT;
END;
IF Done THEN
ReadString (next);
END;
IF STRCMP (next, "up") = 0 THEN
(* read direction of view *)
i := 1;
WHILE Done & (i <= 3) DO
ReadReal (eye.up[i]);
INC (i);
END (* while *);
ELSE
WriteString ("ERROR: viewerdata not correct.");
HALT;
END;
IF Done THEN
ReadString (next);
END;
IF STRCMP (next, "angle") = 0 THEN
(* read viewing angle *)
ReadReal (eye.angle);
ELSE
WriteString ("ERROR: viewing angle not correct.");
HALT;
END;
IF Done THEN
ReadString (next);
END;
IF STRCMP (next, "hither") = 0 THEN
ReadReal (hither);
ELSE
WriteString ("ERROR: viewerdata not correct.");
HALT;
END;
IF Done THEN
ReadString (next);
END;
IF STRCMP (next, "resolution") = 0 THEN
ReadInt (resx);
ReadInt (resy);
ELSE
WriteString ("ERROR: viewerdata not correct.");
HALT;
END;
IF (resx <> Size) OR (resy <> Size) THEN
WriteString ("Using default imagesize for resolution.");
WriteLn;
END;
RETURN;
END get_viewer;
PROCEDURE get_light ();
(***************************************************************************)
(* *)
(* calling procedure: *)
(* createscene *)
(* *)
(* global variable: *)
(* ll *)
(* *)
(* function: *)
(* read data of one light. *)
(* *)
(***************************************************************************)
SCALAR
i: INTEGER; (* counter *)
l: string;
BEGIN
(* increase number of lights and check, if max number is not overcrossed *)
INC (lightcount);
IF (lightcount > MaxLights) THEN
WriteString ("ERROR: too much lights, maximal");
WriteInt (MaxLights,2);
HALT;
END (* if *);
(* read position *)
i := 1;
WHILE Done & (i <= 3) DO
ReadReal (ll[lightcount].pos[i]);
INC (i);
END (* while *);
(* read intensity of light *)
i := 1;
WHILE Done & (i <= 3) DO
ReadReal (ll[lightcount].intensity[i]);
INC (i);
END (* while *);
IF NOT Done THEN
WriteString ("ERROR: lightdata not correct.");
HALT;
END;
END get_light;
PROCEDURE get_background ();
(***************************************************************************)
(* *)
(* calling procedure: *)
(* createscene *)
(* *)
(* global variable: *)
(* back *)
(* *)
(* function: *)
(* read background data. *)
(* *)
(***************************************************************************)
SCALAR
i: INTEGER; (* counter *)
BEGIN
IF back THEN
WriteString ("WARNING: more than one background-color specified.");
WriteLn;
WriteString (" use last one.");
WriteLn;
ELSE
back := TRUE;
END;
i := 1;
WHILE Done & (i <= 3) DO
ReadReal (background[i]);
INC (i);
END (* while *);
IF NOT Done THEN
WriteString ("ERROR: background-data not correct.");
HALT;
END;
END get_background;
PROCEDURE get_material ();
(***************************************************************************)
(* *)
(* calling procedure: *)
(* createscene *)
(* *)
(* global variable: *)
(* *)
(* function: *)
(* read surface-data for the following object. *)
(* *)
(***************************************************************************)
BEGIN
material := TRUE;
IF Done THEN
ReadReal (mat.color[1]);
ReadReal (mat.color[2]);
ReadReal (mat.color[3]);
ReadReal (mat.kdr);
ReadReal (mat.ksr);
ReadReal (mat.shine);
ReadReal (mat.kst);
ReadReal (mat.eta);
END;
IF NOT Done THEN
WriteString ("WARNING: surface-description not correct.");
END;
END get_material;
PROCEDURE get_sphere ();
(***************************************************************************)
(* *)
(* calling procedure: *)
(* createscene *)
(* *)
(* global variable: *)
(* sl *)
(* *)
(* function: *)
(* read sphere-data. *)
(* *)
(***************************************************************************)
SCALAR
i: INTEGER; (* counter *)
BEGIN
(* increase number of spheres and check, if max number is not overcrossed *)
INC (sphcount);
IF (sphcount > MaxSph) THEN
WriteString ("ERROR: too much spheres, maximal");
WriteInt (MaxSph,2);
HALT;
END (* if *);
(* read sphere-data. *)
i := 1;
WHILE (Done) AND (i < 4) DO
ReadReal (sl[sphcount].center[i]);
i := i + 1;
END (* while *);
ReadReal (sl[sphcount].radius);
IF (NOT Done) THEN
WriteString ("ERROR: sphere-format not correct.");
HALT;
END (* if *);
IF NOT material THEN
WriteString ("WARNING: no surface-description for object.");
END;
sl[sphcount].color := mat.color;
sl[sphcount].kdr := mat.kdr;
sl[sphcount].ksr := mat.ksr;
sl[sphcount].kst := mat.kst;
sl[sphcount].shine := mat.shine;
sl[sphcount].eta := mat.eta;
END get_sphere;
PROCEDURE get_poly ();
(***************************************************************************)
(* *)
(* calling procedure: *)
(* createscene *)
(* *)
(* global variable: *)
(* pl *)
(* *)
(* function: *)
(* read polygone-data *)
(* *)
(***************************************************************************)
SCALAR
i: INTEGER; (* counter *)
BEGIN
(* increase number of polygones and check, if max number is not overcrossed*)
INC (polycount);
IF (polycount > MaxPoly) THEN
WriteString ("ERROR: too much polygones, maximal");
WriteInt (MaxPoly,2);
HALT;
END (* if *);
(* read corners of polygones. *)
ReadInt (pl[polycount].vcount);
i := 1;
WHILE (i <= pl[polycount].vcount) & Done DO
ReadReal (pl[polycount].vertices[i][1]);
ReadReal (pl[polycount].vertices[i][2]);
ReadReal (pl[polycount].vertices[i][3]);
INC (i);
END (* while *);
IF NOT Done THEN
WriteString ("ERROR: polygone-data not correct.");
HALT;
END;
IF NOT material THEN
WriteString ("WARNING: no surface-description for object.");
END;
pl[polycount].color := mat.color;
pl[polycount].kdr := mat.kdr;
pl[polycount].ksr := mat.ksr;
pl[polycount].shine := mat.shine;
pl[polycount].kst := mat.kst;
pl[polycount].eta := mat.eta;
END get_poly;
PROCEDURE get_patch ();
(***************************************************************************)
(* *)
(* calling procedure: *)
(* createscene *)
(* *)
(* global variable: *)
(* pl *)
(* *)
(* function: *)
(* read polygone-data. *)
(* *)
(***************************************************************************)
SCALAR
i: INTEGER; (* counter *)
h: REAL;
BEGIN
(* increase number of polygones and check, if max number is not overcrossed*)
INC (polycount);
IF (polycount > MaxPoly) THEN
WriteString ("ERROR: too much polygones, maximal");
WriteInt (MaxPoly,2);
HALT;
END (* if *);
(* read corners of polygones. *)
ReadInt (pl[polycount].vcount);
i := 1;
WHILE (i <= pl[polycount].vcount) & Done DO
ReadReal (pl[polycount].vertices[i][1]);
ReadReal (pl[polycount].vertices[i][2]);
ReadReal (pl[polycount].vertices[i][3]);
(* Normalen ueberlesen *)
ReadReal (h);
ReadReal (h);
ReadReal (h);
INC (i);
END (* while *);
IF NOT Done THEN
WriteString ("ERROR: polygone-data not correct.");
HALT;
END;
IF NOT material THEN
WriteString ("WARNING: no surface-description for object.");
END;
pl[polycount].color := mat.color;
pl[polycount].kdr := mat.kdr;
pl[polycount].ksr := mat.ksr;
pl[polycount].shine := mat.shine;
pl[polycount].kst := mat.kst;
pl[polycount].eta := mat.eta;
END get_patch;
PROCEDURE Read_comment ();
(***************************************************************************)
(* *)
(* calling procedure: *)
(* createscene *)
(* *)
(* function: *)
(* read line untill eoln *)
(* *)
(***************************************************************************)
SCALAR
c: CHAR;
BEGIN
REPEAT
Read (c);
UNTIL (NOT Done) OR (c = EOL);
END Read_comment;
BEGIN (* createscene *)
lightcount := 0;
sphcount := 0;
polycount := 0;
view := FALSE;
back := FALSE;
material := FALSE;
OpenInput (objfile);
IF (Done) THEN
ReadString (t);
ELSE
WriteString ("ERROR: can't open inputfile");
HALT;
END (* if *);
WHILE (Done) DO
IF STRCMP (t, "v") = 0 THEN
get_viewer;
Read_comment;
ELSIF STRCMP (t, "l") = 0 THEN
get_light;
Read_comment;
ELSIF STRCMP (t, "f") = 0 THEN
get_material;
Read_comment;
ELSIF STRCMP (t, "b") = 0 THEN
get_background;
Read_comment;
ELSIF STRCMP (t, "s") = 0 THEN
get_sphere;
Read_comment;
ELSIF STRCMP (t, "p") = 0 THEN
get_poly;
Read_comment ;
ELSIF STRCMP (t, "pp") = 0 THEN
get_patch;
Read_comment ;
ELSIF STRCMP (t, "#") = 0 THEN
Read_comment;
ELSE
WriteString ("ERROR: wrong tag in inputfile.");
HALT;
END;
IF Done THEN
ReadString (t);
END;
END (* while *);
IF NOT back THEN
init_s (background);
END;
IF NOT view THEN
WriteString ("ERROR: no viewer specified.");
HALT;
END;
amb[1] := .035;
amb[2] := .035;
amb[3] := .035;
END createscene;
PROCEDURE planes ();
(***************************************************************************)
(* *)
(* global variable: *)
(* pl, polycount *)
(* *)
(* function: *)
(* calculates out of the first three corners the plain-coefficient of *)
(* the polygone: Ax + By + Cz + D = 0. *)
(* A, B, C represent the normalvector, which is scaled on length 1 *)
(* and adjusting D corresponding *)
(* *)
(***************************************************************************)
SCALAR
x :INTEGER;
VECTOR
v1, v2: Vec; (* edge-vectors *)
i: INTEGER; (* counter *)
r: REAL;
dot: REAL;
BEGIN
x := MaxPoly;
LOAD [1],[1..polycount] (Vecpoly, pl,x);
PARALLEL [1],[1..polycount]
(* calculating coefficients *)
v1[1] := Vecpoly.vertices[2][1] - Vecpoly.vertices[1][1];
v1[2] := Vecpoly.vertices[2][2] - Vecpoly.vertices[1][2];
v1[3] := Vecpoly.vertices[2][3] - Vecpoly.vertices[1][3];
v2[1] := Vecpoly.vertices[3][1] - Vecpoly.vertices[2][1];
v2[2] := Vecpoly.vertices[3][2] - Vecpoly.vertices[2][2];
v2[3] := Vecpoly.vertices[3][3] - Vecpoly.vertices[2][3];
Vecpoly.normal[1] := v1[2] * v2[3] - v1[3] * v2[2];
Vecpoly.normal[2] := v2[1] * v1[3] - v2[3] * v1[1];
Vecpoly.normal[3] := v1[1] * v2[2] - v1[2] * v2[1];
Vecpoly.d := - Vecpoly.normal[1] * Vecpoly.vertices [1][1] -
Vecpoly.normal[2] * Vecpoly.vertices [1][2] -
Vecpoly.normal[3] * Vecpoly.vertices [1][3];
(* check, if points represents a plain *)
IF (Vecpoly.normal[1] = 0.0) AND (Vecpoly.normal[2] = 0.0) AND
(Vecpoly.normal[3] = 0.0) AND (Vecpoly.d = 0.0) THEN
WriteString ("ERROR: polygonpoints represent no plain");
HALT;
END (* if *);
i := 4;
WHILE i <= Vecpoly.vcount DO
IF (ABS(Vecpoly.normal[1]*Vecpoly.vertices[i][1] +
Vecpoly.normal[2]*Vecpoly.vertices[i][2] +
Vecpoly.normal[3]*Vecpoly.vertices[i][3] +
Vecpoly.d) > .001) THEN
WriteString ("ERROR: points represent no plain");
HALT;
END (* if *);
INC (i);
END (* while *);
(* scaling normalvector and adjusting d corresponding. *)
Vecdot_vv (Vecpoly.normal, Vecpoly.normal, dot);
Veccompr_v (1./Sqrt (dot), Vecpoly.normal);
Vecpoly.d := Vecpoly.d / Sqrt(dot);
(* helpvariables for the polygoninintersectionroutine *)
IF (ABS (Vecpoly.normal[3]) <> 0.) THEN
Vecpoly.xy := TRUE;
Vecpoly.xz := FALSE;
Vecpoly.side :=
TRUNC (Vecpoly.normal[3] / ABS (Vecpoly.normal[3]));
ELSIF (ABS (Vecpoly.normal[2]) <> 0.) THEN
Vecpoly.xz := TRUE;
Vecpoly.xy := FALSE;
Vecpoly.side :=
TRUNC (Vecpoly.normal[2] / ABS (Vecpoly.normal[2]));
ELSE
Vecpoly.xy := FALSE;
Vecpoly.xz := FALSE;
Vecpoly.side :=
TRUNC (Vecpoly.normal[1] / ABS (Vecpoly.normal[1]));
END;
ENDPARALLEL;
x := MaxPoly;
STORE [1],[1..polycount] (Vecpoly, pl,x);
END planes;
(*********************** GENERAL HELP-PROCEDURES ***********************)
PROCEDURE initrays ();
(***************************************************************************)
(* *)
(* global variables: *)
(* iray, eye *)
(* *)
(* function: *)
(* Calculation of the initial rays per pixel. *)
(* It would be attempt to distribute as much as possible PEs equal over *)
(* the viewing angle put up by the viewer. As virtual picture-plain there *)
(* would be used the plain in the distance 'eye.at[3]'. This plain is *)
(* only necessary for the calculation of the PE-density per unit. *)
(* The final initial rays depends only on the viewing angle. *)
(* *)
(***************************************************************************)
SCALAR
v: Vec; (* direction *)
up: Vec; (* vertical direction *)
right: Vec; (* horizontal direction *)
VECTOR
t: REAL; (* helping variables *)
BEGIN
Veccomb_ss (1., -1., eye.at, eye.pos, v);
(* calculation of the horizontal rellocation-vector *)
Vecunit_s (v);
vcross_ss (v, eye.up, right);
Vecunit_s (right);
(* calculation of the horizontal rellocation-vector *)
vcross_ss (right,v, up);
Vecunit_s (up);
(* calculation of the direction of initial rays *)
PARALLEL
iray.p := eye.pos;
t := Tan (eye.angle / 360. * PI) *
(-1. + 2. * FLOAT (DIM2 - 1) / FLOAT (Size));
Veccomb_vv (t, 1., right, v, iray.d);
t := Tan (eye.angle / 360. * PI) *
(1. - 2. * FLOAT (DIM1 - 1) / FLOAT (Size));
Veccomb_vv (1., t, iray.d, up, iray.d);
Vecunit_v (iray.d);
ENDPARALLEL;
END initrays;
PROCEDURE round (SCALAR x: REAL): SCALAR INTEGER;
(***************************************************************************)
(* *)
(* function: *)
(* returns rounded integer value *)
(* *)
(***************************************************************************)
SCALAR
z: REAL; (* helping variables *)
BEGIN
z := x * 10. - FLOAT (TRUNC (x)*10);
IF z < 0.0 THEN
IF (z <= -5.) THEN
RETURN (TRUNC (x-1.));
ELSE
RETURN (TRUNC (x));
END (* if *);
ELSE
IF (z >= 5.) THEN
RETURN (TRUNC (x+1.));
ELSE
RETURN (TRUNC (x));
END (* if *);
END (* if *);
END round;
(************************** VECTOR-PROCEDURES *******************************)
PROCEDURE Veccomb_vv (VECTOR a,b: REAL; VECTOR va,vb: Vec; VECTOR VAR vc: Vec);
(***************************************************************************)
(* *)
(* function: *)
(* vectorial vectoraddition: vc = a*va + b*vb *)
(* *)
(* *)
(***************************************************************************)
BEGIN
vc[1] := a * va[1] + b * vb[1];
vc[2] := a * va[2] + b * vb[2];
vc[3] := a * va[3] + b * vb[3];
END Veccomb_vv;
PROCEDURE Veccomb_vs (SCALAR a,b: REAL; VECTOR va: Vec;
SCALAR vb: Vec; VECTOR VAR vc: Vec);
(***************************************************************************)
(* *)
(* function: *)
(* vectoraddition: vc = a*va + b*vb *)
(* *)
(* *)
(***************************************************************************)
BEGIN
vc[1] := a * va[1] + b * vb[1];
vc[2] := a * va[2] + b * vb[2];
vc[3] := a * va[3] + b * vb[3];
END Veccomb_vs;
PROCEDURE Veccomb_ss (SCALAR a,b: REAL; SCALAR va: Vec;
SCALAR vb: Vec; SCALAR VAR vc: Vec);
(***************************************************************************)
(* *)
(* function: *)
(* vectoraddition: vc = a*va + b*vb *)
(* *)
(* *)
(***************************************************************************)
BEGIN
vc[1] := a * va[1] + b * vb[1];
vc[2] := a * va[2] + b * vb[2];
vc[3] := a * va[3] + b * vb[3];
END Veccomb_ss;
PROCEDURE Vecdot_vv (VECTOR va, vb: Vec; VECTOR VAR d: REAL);
(***************************************************************************)
(* *)
(* function: *)
(* vectorial vectormultiplication: vc = a*va + vb *)
(* *)
(***************************************************************************)
BEGIN
d := va[1]*vb[1] + va[2]*vb[2] + va[3]*vb[3];
END Vecdot_vv;
PROCEDURE Vecdot_sv (SCALAR va: Vec; VECTOR vb: Vec; VECTOR VAR d: REAL);
(***************************************************************************)
(* *)
(* function: *)
(* vectorial vectormultiplication: vc = a*va + vb *)
(* *)
(***************************************************************************)
BEGIN
d := va[1]*vb[1] + va[2]*vb[2] + va[3]*vb[3];
END Vecdot_sv;
PROCEDURE Vecdot_ss (SCALAR va, vb: Vec; SCALAR VAR d: REAL);
(***************************************************************************)
(* *)
(* function: *)
(* scalar vectormultiplication: vc = a*va + vb *)
(* *)
(***************************************************************************)
BEGIN
d := va[1]*vb[1] + va[2]*vb[2] + va[3]*vb[3];
END Vecdot_ss;
PROCEDURE vcross_ss (SCALAR va, vb: Vec; SCALAR VAR vc: Vec);
(***************************************************************************)
(* *)
(* function: *)
(* calculation of the crossproduct of va and vb. *)
(* *)
(***************************************************************************)
BEGIN
vc[1] := va[2] * vb[3] - va[3] * vb[2];
vc[2] := va[3] * vb[1] - va[1] * vb[3];
vc[3] := va[1] * vb[2] - va[2] * vb[1];
END vcross_ss;
PROCEDURE Veccompr_v (VECTOR a: REAL; VECTOR VAR va: Vec);
(***************************************************************************)
(* *)
(* function: *)
(* stretching the vector va: va = a*va *)
(* *)
(***************************************************************************)
BEGIN
va[1] := a * va[1];
va[2] := a * va[2];
va[3] := a * va[3];
END Veccompr_v;
PROCEDURE Veccompr_s (SCALAR a: REAL; SCALAR VAR va: Vec);
(***************************************************************************)
(* *)
(* function: *)
(* stretching the vector va: va = a*va *)
(* *)
(***************************************************************************)
BEGIN
va[1] := a * va[1];
va[2] := a * va[2];
va[3] := a * va[3];
END Veccompr_s;
PROCEDURE Vecunit_v (VECTOR VAR va: Vec);
(***************************************************************************)
(* *)
(* function: *)
(* vectorial calculation of the unitvector of va *)
(* *)
(***************************************************************************)
VECTOR
d: REAL;
BEGIN
Vecdot_vv (va, va, d);
d := 1./Sqrt(d);
Veccompr_v (d, va);
END Vecunit_v;
PROCEDURE Vecunit_s (SCALAR VAR va: Vec);
(***************************************************************************)
(* *)
(* function: *)
(* vectorial calculation of the unitvector of va *)
(* *)
(***************************************************************************)
SCALAR
d: REAL;
BEGIN
Vecdot_ss (va, va, d);
d := 1./Sqrt(d);
Veccompr_s (d, va);
END Vecunit_s;
PROCEDURE dist (SCALAR va: Vec; VECTOR vb: Vec; VECTOR VAR d: REAL);
(***************************************************************************)
(* *)
(* function: *)
(* calculation of the distance between two points *)
(* *)
(***************************************************************************)
BEGIN
d:= Sqrt ((va[1]-vb[1])**2 + (va[2]-vb[2])**2 + (va[3]-vb[3])**2);
END dist;
(************************* COLOR-PROcEDUREs **********************************)
PROCEDURE init_s (SCALAR VAR col: RGB);
(***************************************************************************)
(* *)
(* function: *)
(* col get the color black *)
(* *)
(***************************************************************************)
BEGIN
col[1] := 0.;
col[2] := 0.;
col[3] := 0.;
END init_s;
PROCEDURE init_v (VECTOR VAR col: RGB);
(***************************************************************************)
(* *)
(* function: *)
(* vector col get the color black *)
(* *)
(***************************************************************************)
BEGIN
col[1] := 0.;
col[2] := 0.;
col[3] := 0.;
END init_v;
PROCEDURE coladd_vv (VECTOR a: REAL; VECTOR va, vb: RGB; VECTOR VAR vc: RGB);
(***************************************************************************)
(* *)
(* function: *)
(* addition of colorvectors. *)
(* *)
(***************************************************************************)
BEGIN
vc[1] := a * va[1] + vb[1];
vc[2] := a * va[2] + vb[2];
vc[3] := a * va[3] + vb[3];
END coladd_vv;
PROCEDURE coladd_sv (VECTOR a: REAL; SCALAR va: RGB; VECTOR vb: RGB;
VECTOR VAR vc: RGB);
(***************************************************************************)
(* *)
(* function: *)
(* addition of colorvectors. *)
(* *)
(***************************************************************************)
BEGIN
vc[1] := a * va[1] + vb[1];
vc[2] := a * va[2] + vb[2];
vc[3] := a * va[3] + vb[3];
END coladd_sv;
PROCEDURE colmult_sv (SCALAR va: RGB; VECTOR vb: RGB; VECTOR VAR vc: RGB);
(***************************************************************************)
(* *)
(* function: *)
(* multiplication of colorvectors. *)
(* *)
(***************************************************************************)
BEGIN
vc[1] := va[1] * vb[1];
vc[2] := va[2] * vb[2];
vc[3] := va[3] * vb[3];
END colmult_sv;
(*************************** OBJECT-PROCEDURES ******************************)
PROCEDURE sphere_normal (VECTOR s: Sphere; VECTOR p: Vec; VECTOR VAR n: Vec);
(***************************************************************************)
(* *)
(* function: *)
(* calculating the normalvector of sphere s at point p. *)
(* *)
(***************************************************************************)
BEGIN
Veccomb_vv (-1.,1., s.center, p, n);
Vecunit_v (n);
END sphere_normal;
PROCEDURE poly_normal (VECTOR poly: Polygon; VECTOR VAR n: Vec);
(***************************************************************************)
(* *)
(* function: *)
(* calculating the normalvector of polygon poly *)
(* *)
(***************************************************************************)
BEGIN
n[1] := poly.normal[1];
n[2] := poly.normal[2];
n[3] := poly.normal[3];
END poly_normal;
PROCEDURE sphere_intersect (VECTOR p, d: Vec; VECTOR VAR t: REAL;
VECTOR VAR typ: ObjTyp; VECTOR VAR s: Sphere;
VECTOR VAR hit: BOOLEAN);
(***************************************************************************)
(* *)
(* global variable: *)
(* sl *)
(* *)
(* function: *)
(* Calculates the nearest intersection between ray with origin 'p' and *)
(* direction 'd'. The distance between intersection and origin has to be *)
(* less than 't'. *)
(* The whole spherelist 'sl' would be passed and each new intersection *)
(* would be checked, if it lays before the other intersections. At the end*)
(* get s the cutting sphere and t the multiplicator of the directionvector*)
(* of which the intersection can be calculated. *)
(* *)
(***************************************************************************)
SCALAR
i: INTEGER; (* counter *)
VECTOR
v: Vec; (* vector to center of sphere *)
b, vdot, disc, t1, t2: REAL;
inters: BOOLEAN;
BEGIN
FOR i := 1 TO sphcount DO
Veccomb_vs (-1., 1., p, sl[i].center, v);
Vecdot_vv (v, d, b);
Vecdot_vv (v, v, vdot);
disc := b*b - vdot + sl[i].radius**2;
IF (disc <= 0.) THEN (* no cut *)
inters := FALSE;
ELSE
disc := Sqrt (disc);
t2 := b + disc;
IF (t2 < RayEps) THEN
inters := FALSE;
ELSE
inters := TRUE;
t1 := b - disc;
IF (t1 < RayEps) THEN
t1 := t2;
END (* if *);
END (* if *)
END (* if *);
IF (inters) & (t > t1) THEN (* new intersection lays before *)
(* the old one *)
t := t1;
typ := sph;
s := sl[i];
END (* if *);
END (* for *);
IF typ = sph THEN
hit := TRUE;
ELSE
hit := FALSE;
END (* if *);
END sphere_intersect;
PROCEDURE poly_intersect (VECTOR p, d: Vec; VECTOR VAR tmax: REAL;
VECTOR VAR typ: ObjTyp; VECTOR VAR pn: Polygon;
VECTOR VAR hit: BOOLEAN);
(***************************************************************************)
(* *)
(* subprocedure: *)
(* in_out_test *)
(* *)
(* global variable: *)
(* pl *)
(* *)
(* function: *)
(* calculates intersection of ray, which run from p in direction d. *)
(* *)
(***************************************************************************)
SCALAR
i: INTEGER;
VECTOR
cosdn: REAL; (* pointproduct of normalvector *)
(* and falling-in ray *)
cospn: REAL; (* pointproduct of normalvector *)
(* and origin *)
t: REAL; (* multiplicator of direction- *)
(* vector *)
point: Vec; (* meeting-point *)
inside: BOOLEAN; (* TRUE: point lay inside polygon*)
(* FALSE: point lay outside *)
PROCEDURE in_out_test (VECTOR p: Vec; SCALAR plentry: INTEGER;
VECTOR VAR inpoly: BOOLEAN);
(***************************************************************************)
(* *)
(* calling procedure: *)
(* poly_intersect *)
(* *)
(* global variable: *)
(* pl *)
(* *)
(* function: *)
(* Check, if point 'p' lay inside polygone at position 'plentry'. *)
(* Send one ray in one direction and count the intersections with *)
(* polygone-edges. If the number is even, then the point lay inside, else *)
(* outside. *)
(* *)
(***************************************************************************)
SCALAR
i: INTEGER; (* counter *)
v1, v2: Vec; (* current corner-points of *)
BEGIN
i := 1;
inpoly := TRUE;
WHILE inpoly & (i <= pl[plentry].vcount) DO
IF i = pl[plentry].vcount THEN
v2 := pl[plentry].vertices[1];
ELSE
v2 := pl[plentry].vertices[i+1];
END (* if *);
v1 := pl[plentry].vertices[i];
IF pl[plentry].xy THEN
(* Projektion auf die x,y - Ebene *)
IF (((v1[2] - v2[2]) * p[1] + (v2[1] -v1[1]) * p[2] +
(v1[1]*v2[2]-v2[1]*v1[2]))*FLOAT (pl[plentry].side)<0.) THEN
inpoly := FALSE;
END (* if *);
ELSIF pl [plentry].xz THEN
(* projection to x,z - plain *)
IF (((v1[3] - v2[3]) * p[1] + (v2[1] -v1[1]) * p[3] +
(v1[1]*v2[3]-v2[1]*v1[3]))*FLOAT (pl[plentry].side)<0.) THEN
inpoly := FALSE;
END (* if *);
ELSE
(* projection to y,z - plain *)
IF (((v1[3] - v2[3]) * p[2] + (v2[2] -v1[2]) * p[3] +
(v1[2]*v2[3]-v2[2]*v1[3]))*FLOAT (pl[plentry].side)<0.) THEN
inpoly := FALSE;
END (* if *);
END;
INC (i);
END;
END in_out_test;
BEGIN
FOR i := 1 TO polycount DO
Vecdot_sv (pl[i].normal, d, cosdn);
IF cosdn <> 0. THEN (* ray not parallel with plain *)
Vecdot_sv (pl[i].normal, p, cospn);
t := (-cospn - pl[i].d) / cosdn;
IF (t > RayEps) THEN (* hit lay before ray-origin *)
Veccomb_vv (1., t, p, d, point);
in_out_test (point, i, inside);
IF inside & (t < tmax) THEN
tmax := t;
pn := pl[i];
typ := poly
END;
END;
END;
END;
IF typ = poly THEN
hit := TRUE;
END;
END poly_intersect;
(********************* ACTUAL RAYTRACE PROCEDURES *********************)
PROCEDURE trace (SCALAR depth: INTEGER; VECTOR VAR color: RGB;
VECTOR VAR r: Ray; VECTOR weight: REAL);
(***************************************************************************)
(* *)
(* global variable: *)
(* color *)
(* *)
(* function: *)
(* Calculates intersection between 'r' and the nearest object. *)
(* Call of shadowing routine. *)
(* *)
(***************************************************************************)
VECTOR
t: REAL; (* parameter for calculating *)
(* intersection *)
obj: Object;
hit: BOOLEAN;
hitpoint: Vec; (* hitpoint *)
normal: Vec; (* normalvec at hitpoint *)
newcolor: RGB; (* color at hitpoint *)
BEGIN
t := 10000.;
obj.typ := none;
IF sphcount > 0 THEN
sphere_intersect (r.p, r.d, t, obj.typ,
obj.sph, hit);
END;
IF polycount > 0 THEN
poly_intersect (r.p, r.d, t, obj.typ,
obj.poly, hit);
END;
IF (obj.typ = sph) OR (obj.typ = poly) THEN
(* all rays which hit object *)
Veccomb_vv (1., t, r.p, r.d, hitpoint);
IF obj.typ = sph THEN
sphere_normal (obj.sph, hitpoint, normal);
ELSE
poly_normal (obj.poly, normal);
END;
shade (depth, obj, hitpoint, r.d, normal, color, weight);
ELSE
color := background;
END (* if *);
END trace;
PROCEDURE shade (SCALAR depth: INTEGER; VECTOR obj: Object;
VECTOR point, d, normal: Vec;
VECTOR VAR color: RGB; VECTOR weight: REAL);
(***************************************************************************)
(* *)
(* global variable: *)
(* ll, sl, pl, amb, color *)
(* *)
(* function: *)
(* Calculation of viewing ray d, which hit object in obj.sph respectively *)
(* obj.poly (look obj.typ) at point p. Phong-Shading is used as shadowing *)
(* model. *)
(* *)
(***************************************************************************)
SCALAR
i: INTEGER; (* counter *)
VECTOR
lcolor: RGB; (* color by light *)
l: Vec; (* shading-ray to light *)
dst: REAL; (* distance of light *)
sh: BOOLEAN; (* TRUE: point is lighted *)
(* FALSE: object lay between *)
(* point and light *)
cosln: REAL; (* pointproduct of shading-ray *)
(* and normalvec *)
cosnh: REAL; (* pointproduct of halfway-vec *)
(* and normalvec *)
h: Vec; (* halfway-vec between shading *)
(* and falling-in ray *)
hdot: REAL; (* pointproduct of halway-vec *)
f: RGB; (* calculated Fresnel-values *)
oldweight: REAL; (* weighting factor of falling-in*)
(* ray *)
spec: BOOLEAN;
trans: BOOLEAN; (* TRUE: ray would be transm. *)
(* FALSE: ray would be total refl*)
newray: Ray;
specweight: REAL;
transweight: REAL;
newweight: REAL;
ksr, kst, mult: REAL;
newcolor: RGB;
BEGIN
(* colorquotas by ambient light and light-origins can be calculated direct*)
(* for each ray. *)
(* Phong - Shading: colorquota := amb. + diffus and spec. refl. light *)
IF obj.typ = sph THEN
colmult_sv (amb, obj.sph.color, color);
ELSE
colmult_sv (amb, obj.poly.color, color);
END;
FOR i := 1 TO lightcount DO
init_v (lcolor);
Veccomb_vs (-1., 1., point, ll[i].pos, l);
dist (ll[i].pos, point, dst);
Vecunit_v (l);
Vecdot_vv ( l, normal, cosln);
shadow (point, l, dst, sh);
IF sh & (cosln > 0.) THEN (* no object in this direction *)
(* halfway-vec between normalvec and directionvec for spec. reflection *)
Veccomb_vv (-1., 1., d, l, h);
Vecdot_vv (h, h, hdot);
Veccompr_v (1. / Sqrt(hdot), h);
Vecdot_vv (normal, h, cosnh);
IF obj.typ = sph THEN
(* diffuser quota *)
colmult_sv (ll[i].intensity,obj.sph.color,lcolor);
coladd_vv (obj.sph.kdr*cosln,lcolor, color,color);
(* spec. quota *)
coladd_sv
(obj.sph.ksr*cosnh**obj.sph.shine, ll[i].intensity,
color, color);
ELSE
(* diffuse quota *)
colmult_sv (ll[i].intensity,obj.poly.color,lcolor);
coladd_vv (obj.poly.kdr*cosln,lcolor, color,color);
(* spec. quota *)
coladd_sv
(obj.poly.ksr*cosnh**obj.poly.shine, ll[i].intensity,
color, color);
END;
END;
END (* for *);
(* max. depth still not achieved --> follow spec. and trans. ray further *)
IF (depth < MaxDepth) THEN
IF obj.typ = sph THEN
specweight := weight * obj.sph.ksr;
transweight := weight * obj.sph.kst;
ksr := obj.sph.ksr;
kst := obj.sph.kst;
ELSE
specweight := weight * obj.poly.ksr;
transweight := weight * obj.poly.kst;
ksr := obj.poly.ksr;
kst := obj.poly.kst;
END;
newray.p := point;
init_v (newcolor);
spec := FALSE;
trans := FALSE;
IF (specweight > MinWeight) THEN
(* weight of spec. refl. ray relevant --> calculate spec. refl. rays *)
specdir (d, normal, newray.d);
newweight := specweight;
spec := TRUE;
mult := ksr;
ELSIF (transweight > MinWeight) THEN
transdir (obj, d, normal, newray.d, trans);
IF NOT trans THEN
specdir (d, normal, newray.d);
trans := TRUE;
END;
newweight := transweight;
mult := kst;
END;
IF spec OR trans THEN
trace (depth + 1, newcolor, newray, newweight);
coladd_vv (mult, newcolor, color, color);
END;
IF spec AND (transweight > MinWeight) THEN
(* weight of transm. ray relevant --> calculate transm. rays *)
transdir (obj, d, normal, newray.d, trans);
IF NOT trans THEN
specdir (d, normal, newray.d);
END;
init_v (newcolor);
trace (depth + 1, newcolor, newray, newweight);
coladd_vv (kst, newcolor, color,color);
END (* if *);
END (* if *);
END shade;
PROCEDURE shadow (VECTOR p, d: Vec; VECTOR tmax: REAL; VECTOR VAR b: BOOLEAN);
(***************************************************************************)
(* *)
(* function: *)
(* Check, if shading-ray running from p hit another object, before he *)
(* reach the light in distance tmax. *)
(* *)
(***************************************************************************)
VECTOR
hit: BOOLEAN;
t: REAL;
hitobj: Object;
BEGIN
t := tmax;
hitobj.typ := none;
hit := FALSE;
IF sphcount > 0 THEN
sphere_intersect (p, d, t, hitobj.typ, hitobj.sph, hit);
END;
IF (NOT hit) OR (t > tmax - RayEps) THEN
IF polycount > 0 THEN
poly_intersect (p, d, t, hitobj.typ, hitobj.poly, hit);
END;
IF (NOT hit) OR (t > tmax - RayEps) THEN
b := TRUE; (* ray reach light *)
ELSE
b := FALSE;
END;
ELSE
b := FALSE; (* ray hit object *)
END (* if *);
END shadow;
PROCEDURE specdir (VECTOR d, n: Vec; VECTOR VAR r: Vec);
(***************************************************************************)
(* *)
(* function: *)
(* Calculates spec. reflect. ray of falling-in ray d at position with *)
(* normalvec n. *)
(* *)
(***************************************************************************)
VECTOR
dot: REAL;
BEGIN
Vecdot_vv (d, n, dot);
Veccomb_vv (-2.*dot, 1., n, d, r);
END specdir;
PROCEDURE transdir (VECTOR obj: Object; VECTOR d, n: Vec;
VECTOR VAR r: Vec; VECTOR VAR b: BOOLEAN);
(***************************************************************************)
(* *)
(* function: *)
(* Calculates Transmission-direction of ray d through obj.sph resp. *)
(* obj.poly at position with normalvec n. *)
(* *)
(***************************************************************************)
VECTOR
eta, c1, cs2: REAL;
BEGIN
IF obj.typ = sph THEN
eta := obj.sph.eta;
ELSE
eta := obj.poly.eta;
END;
Vecdot_vv (d, n, c1);
IF c1 > 0. THEN
eta := 1./eta;
ELSE
c1 := -c1;
END;
cs2 := 1. - eta**2 * (1. - c1**2);
IF (cs2 < 0.) THEN
b := FALSE;
END;
Veccomb_vv (eta, eta*c1-Sqrt(cs2),d, n, r);
b := TRUE;
END transdir;
BEGIN (* main *)
(* planes initrays round Veccomb_vv Veccomb_vs Veccomb_ss
Vecdot_vv Vecdot_sv Vecdot_ss vcross_ss Veccompr_v Veccompr_s
Vecunit_v Vecunit_s dist coladd_vv coladd_sv colmult_sv
sphere_normal poly_normal sphere_intersect
poly_intersect in_out_test trace shade shadow specdir transdir *)
WriteString ("Objectfile (without extension '.nff') : ");
ReadString (inputf);
objfile := strcat( inputf , ".nff");
picfile := strcat( inputf , ".out.ppm");
WriteLn;
WriteString( " Objectfile : " );
WriteString( objfile ); WriteLn;
WriteString( " Outputfile : " );
WriteString( picfile ); WriteLn; WriteLn;
createscene;
IF polycount > 0 THEN
planes;
END;
initrays;
PARALLEL
trace (0, color, iray, 1.);
ENDPARALLEL;
STORE (color, pixfield);
OpenOutput (picfile);
WriteString( "P6"); WriteLn;
WriteInt (Size, 1); WriteLn;
WriteInt (Size, 1); WriteLn;
WriteInt (255, 1) ; WriteLn;
FOR i := 1 TO Size DO
FOR j := 1 TO Size DO
aktu_color := pixfield[i][j];
(* normalize color *)
m := aktu_color[1];
IF aktu_color[2] > m THEN
m := aktu_color[2];
END;
IF aktu_color[3] > m THEN
m := aktu_color[3];
END;
IF m > 1. THEN
aktu_color[1] := aktu_color[1] / m;
aktu_color[2] := aktu_color[2] / m;
aktu_color[3] := aktu_color[3] / m;
END;
Write ( CHR(round (aktu_color[1]*255.)));
Write ( CHR(round (aktu_color[2]*255.)));
Write ( CHR(round (aktu_color[3]*255.)));
END;
END;
CloseOutput;
END sys_trace.